Traveling Salesperson Problem: UPS

Our group identified a vehicle routing problem (VRP) for the shipping and receiving company UPS. Our goal is to identify the most effective routing for UPS delivery drivers within the Dallas locale

Introduction UPS is a leading global package delivery company that operates a vast fleet of vehicles to deliver packages to millions of customers every day. Optimizing routes for this fleet is a critical task that has a significant impact on the company's operations and bottom line. The Vehicle Routing Problem (VRP) is a well-known optimization problem that involves finding the most efficient routes for a fleet of vehicles to deliver goods to a set of customers, considering factors such as distance, delivery time windows, and vehicle capacity constraints.

We used 12 Chipotle locations in the Dallas area to conveniently use the UPS fleet around the distribution center we chose.

Google OR-Tools Google OR-Tools is an open-source software suite for operations research that includes a collection of tools for solving a variety of optimization problems, including the VRP. OR-Tools provides a powerful and flexible platform for solving the VRP, and it has been used by a variety of organizations to improve their transportation operations.

Advantages of using OR-Tools *Flexibility: OR-Tools is a highly flexible platform that can be used to model a wide range of VRP problems, considering a variety of constraints and objective functions.

*Speed: OR-Tools is designed to solve large-scale optimization problems quickly, making it an ideal solution for organizations like UPS that need to optimize routes in real-time.

*Open-source: As an open-source platform, OR-Tools can be freely downloaded and used, which can help reduce costs for organizations like UPS.

*Ease of use: OR-Tools provides a user-friendly interface for solving VRP problems, making it accessible to a wide range of users, including those without extensive technical expertise.g Salesperson Problem: UPS Our group identified a vehicle routing problem (VRP) for the shipping and receiving company UPS. Our goal is to identify the most effective routing for UPS delivery drivers within the Dallas locale

1.Import the modules

In [1]:
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
from sklearn.manifold import MDS

import requests
import json
from urllib.request import urlopen
import matplotlib.pyplot as plt
import math

2.Create Data model
   a. Function defined to take the data which is locations and number of verhicles as parameters.
   b.returns the data dictionary updated with the values. currently the default no of vehicles is 4

In [2]:
def create_data_model(data,num_veh=4):
    data['num_vehicles'] = num_veh
    data['depot'] = 0
    return data

3.Create Data:
  1. Create the data dictionary with the first element as the starting point in this case UPS distribution center
  2. All other values are lat, lon coordinates of the chipotle destination (Stops)

In [3]:
def create_data():
    """Creates the data."""
    data = {}
    data['API_key'] = 'AIzaSyAMDZOhS7STXhHaDyWsoLYiF1nGDIICAh0'
    data['addresses'] = ['32.8846928%2C-96.8817665',
                    '32.90781142%2C-96.76975673',
                    '32.91141547%2C-96.80302352',
                    '32.92135201%2C-96.83873005',
                    '32.94971005%2C-96.76948706',
                    '32.99912051%2C-96.79672475',
                    '32.81639867%2C-96.75316057',
                    '32.78812%2C-96.81001',
                    '32.8273184%2C-96.8455045',
                    '32.79815888%2C-96.80163622',
                    '32.86112926%2C-96.85722228',
                    '32.81777982%2C-96.81203523',
                    '32.81777982%2C-96.81203523'
                      ]
    return data

Plot Using Folium

In [5]:
import folium

# Create a map centered at the given latitude and longitude
map_center = [32.8846928, -96.8817665] # latitude and longitude of London
m = folium.Map(location=map_center, zoom_start=13)

# Add a marker for the given latitude and longitude
#marker_location = [51.5074, -0.1278] # latitude and longitude of the marker
locations =[ [32.8846928, -96.8817665], [32.90781142, -96.76975673], [32.91141547, -96.80302352], [32.92135201, -96.83873005], [32.94971005, -96.76948706], [32.99912051, -96.79672475], [32.81639867, -96.75316057], [32.78812, -96.81001], [32.8273184, -96.8455045], [32.79815888, -96.80163622], [32.86112926, -96.85722228],[32.81777982, -96.81203523],[32.81777982, -96.81203523]]
tooltip_text = ['UPS Distrubution center','A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L'] # tooltip text for the marker
for i in range(len(locations)):    
    folium.Marker(location=locations[i], tooltip=tooltip_text[i]).add_to(m)

# Display the map
m
Out[5]:
Make this Notebook Trusted to load map: File -> Trust Notebook

PLOT 2D locations for visualization of Routing purposes

In [6]:
# Some tooling to ease drawing
height = 7

def plot_location(location, axes, color, location_number):
  axes.scatter(
      location[0],
      location[1],
      s=1000,
      facecolors='white',
      edgecolors=color,
      linewidths=2)
  axes.scatter(
      location[0],
      location[1],
      s=400,
      marker=f'${location_number}$',
      edgecolors=color,
      facecolors=color)

# A diagram of the city is shown below, with the company location marked in black and the locations to visit in blue.
def plot_locations(locations):
  fig, axes = plt.subplots(figsize=(1.7 * height, height))
  axes.grid(True)
  axes.set_xticks(list(set([x for (x, y) in locations])))
  axes.set_xticklabels([])
  axes.set_yticks(list(set([y for (x, y) in locations])))
  axes.set_yticklabels([])
  axes.set_axisbelow(True)
  for (i, location) in enumerate(locations):
    color = 'blue' if i else 'black'
    plot_location(location, axes, color, i)

# Create the 2D Locations
loc =[(4831.7372137462935, 10575.724459700612), (15313.907201992086, 13149.286926238412), (12200.713167728254, 13550.489411833396), (8859.201823223353, 14656.624050979197), (15339.143632534458, 17813.438227975257), (12790.168693410264, 23313.80569412894), (16867.01966952093, 2973.228737513024), (11546.899924655138, 174.74997901984662), (8225.23090732053, 4188.811989140558), (12330.54036427185, 942.777138692521), (7128.650119369224, 7952.633543254063), (6346.614622029878, 15508.561924872702), (11357.373554542019, 3126.9782173980548)]
#print(loc)
plot_locations(loc)
plt.show()

4.Creating Distance Matrix: To compute the distance matrix, we would like to send a single request containing all 13 addresses as both the origin and destination addresses. However, we can't because this would require 13x13=269 origin-destination pairs, while the API is restricted to 100 such pairs per request. So we need to make multiple requests.

Since each row of the matrix contains 13 entries, we can compute at most seven rows per request (requiring 7x13=91 pairs). We can compute the whole matrix in three requests, which return 7 rows, and 6 rows.

The following code computes the distance matrix as follows:

Divide the 13 addresses into two groups of seven addresses and one group of six addresses. For each group, build and send a request for the origin addresses in the group and all destination addresses. See Build and send a request. Use the response to build the corresponding rows of the matrix, and concatenate the rows (which are just Python lists). See Build rows of the distance matrix.

In [7]:
def create_distance_matrix(data):
    addresses = data["addresses"]
    API_key = data["API_key"]
    # Distance Matrix API only accepts 100 elements per request, so get rows in multiple requests.
    max_elements = 100
    num_addresses = len(addresses) # 13 in this example.
    # Maximum number of rows that can be computed per request (7 in this example).
    max_rows = max_elements // num_addresses
    # num_addresses = q * max_rows + r (q = 1 and r = 6 in this example).
    q, r = divmod(num_addresses, max_rows)
    print(f'q,r: {q,r,len(addresses)}')
    dest_addresses = addresses
    distance_matrix = []
    # Send q requests, returning max_rows rows per request.
    for i in range(q):
        origin_addresses = addresses[i * max_rows: (i + 1) * max_rows]
        response = send_request(origin_addresses, dest_addresses, API_key)
        distance_matrix += build_distance_matrix(response)

    # Get the remaining remaining r rows, if necessary.
    if r > 0:
        origin_addresses = addresses[q * max_rows: q * max_rows + r]
        response = send_request(origin_addresses, dest_addresses, API_key)
        distance_matrix += build_distance_matrix(response)
    return distance_matrix

def send_request(origin_addresses, dest_addresses, API_key):
    """ Build and send request for the given origin and destination addresses."""
    def build_address_str(addresses):
        # Build a pipe-separated string of addresses
        address_str = ''
        for i in range(len(addresses) - 1):
          address_str += addresses[i] + '|'
        address_str += addresses[-1]
        return address_str

    request = 'https://maps.googleapis.com/maps/api/distancematrix/json?units=imperial'
    origin_address_str = build_address_str(origin_addresses)
    dest_address_str = build_address_str(dest_addresses)
    request = request + '&origins=' + origin_address_str + '&destinations=' + \
                       dest_address_str + '&key=' + API_key
    jsonResult = urlopen(request).read()
    response = json.loads(jsonResult)
    return response

def build_distance_matrix(response):
    distance_matrix = []
    for row in response['rows']:
        row_list = [row['elements'][j]['distance']['value'] for j in range(len(row['elements']))]
        distance_matrix.append(row_list)
    return distance_matrix

5.Print Solution an PLot it:
  a.This piece of code prints the solution derived from the manager object
  b. The other function plot_solution plots the route in google stype (Blue, red, yellow ) colors

In [8]:
def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f'Objective: {solution.ObjectiveValue()}')
    max_route_distance = 0
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        plan_output = 'Route for vehicle {}:\n'.format(vehicle_id)
        route_distance = 0
        while not routing.IsEnd(index):
            plan_output += ' {} -> '.format(manager.IndexToNode(index))
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(
                previous_index, index, vehicle_id)
        plan_output += '{}\n'.format(manager.IndexToNode(index))
        plan_output += 'Distance of the route: {}m\n'.format(route_distance)
        print(plan_output)
        max_route_distance = max(route_distance, max_route_distance)
    print('Maximum of the route distances: {}m'.format(max_route_distance))
    return max_route_distance

# a diagram of the solution using a Google colorscheme
def plot_solution(locations, manager, routing, solution, loc, marker_size):
  height = 8
  fig, axes = plt.subplots(figsize=(1.7 * height, height))
  axes.grid(True)
  axes.set_xticks(list(set([x for (x, y) in locations])))
  axes.set_xticklabels([])
  axes.set_yticks(list(set([y for (x, y) in locations])))
  axes.set_yticklabels([])
  axes.set_axisbelow(True)
  max_route_distance = 0
  google_colors = [
      r'#4285F4', r'#EA4335', r'#FBBC05', r'#34A853', r'#101010', r'#FFFFFF'
  ]
  for vehicle_id in range(manager.GetNumberOfVehicles()):
    previous_index = routing.Start(vehicle_id)
    while not routing.IsEnd(previous_index):
      index = solution.Value(routing.NextVar(previous_index))
      start_node = manager.IndexToNode(previous_index)
      end_node = manager.IndexToNode(index)
      start = locations[start_node]
      end = locations[end_node]
      delta_x = end[0] - start[0]
      delta_y = end[1] - start[1]
      delta_length = math.sqrt(delta_x**2 + delta_y**2)
      if delta_length == 0.0:
            previous_index = index
            continue
      unit_delta_x = delta_x / delta_length
      unit_delta_y = delta_y / delta_length
      axes.arrow(
          start[0] + (marker_size / 2) * unit_delta_x,
          start[1] + (marker_size / 2) * unit_delta_y,
          (delta_length - marker_size) * unit_delta_x,
          (delta_length - marker_size) * unit_delta_y,
          head_width=20,
          head_length=20,
          facecolor=google_colors[vehicle_id],
          edgecolor=google_colors[vehicle_id],
          length_includes_head=True,
          width=5)
      previous_index = index
      node_color = 'black' if routing.IsEnd(
          previous_index) else google_colors[vehicle_id]
      plot_location(end, axes, node_color, end_node)

MainProgram
  the main program fucntion with OR tools calculates the distance dimentions and returns the solution

In [9]:
def opt_vrp(num_veh):
    """Entry point of the program"""
    # Create the data.
    data = create_data()
    addresses = data['addresses']
    API_key = data['API_key']
    data['distance_matrix']= create_distance_matrix(data)
    #print(data['distance_matrix'])
    # Instantiate the data problem.
    create_data_model(data,num_veh)
    #print(data)

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),
                                           data['num_vehicles'], data['depot'])
    print(manager)
    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback.
    def distance_callback(from_index, to_index):
        """Returns the distance between the two nodes."""
        # Convert from routing variable Index to distance matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data['distance_matrix'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Distance constraint.
    dimension_name = 'Distance'
    routing.AddDimension(
        transit_callback_index,
        0,  # no slack
        300000,  # vehicle maximum travel distance
        True,  # start cumul to zero
        dimension_name)
    distance_dimension = routing.GetDimensionOrDie(dimension_name)
    distance_dimension.SetGlobalSpanCostCoefficient(100)

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)
    #print(f"solution:{solution}")

    # Print solution on console.
    if solution:
        maxroutedistance = print_solution(data, manager, routing, solution)
        plot_solution(loc, manager, routing, solution, loc,  50)
        plt.show()
        return maxroutedistance
    else:
        print('No solution found !')
        
    
In [10]:
solution = opt_vrp(4)
q,r: (1, 6, 13)
<ortools.constraint_solver.pywrapcp.RoutingIndexManager; proxy of <Swig Object of type 'operations_research::RoutingIndexManager *' at 0x7f16e59b3180> >
Objective: 4240090
Route for vehicle 0:
 0 ->  3 ->  4 ->  1 -> 0
Distance of the route: 35110m

Route for vehicle 1:
 0 ->  8 ->  7 ->  9 ->  12 -> 0
Distance of the route: 32494m

Route for vehicle 2:
 0 ->  2 ->  5 -> 0
Distance of the route: 40940m

Route for vehicle 3:
 0 ->  6 ->  11 ->  10 -> 0
Distance of the route: 37546m

Maximum of the route distances: 40940m

Sensitivity Analysis   a. Ran a sensitvity analysis for numberofvehicles and maxdistance. As earlier the API we have kept the number of vehicles as constant and find out the optimzed route.   b. the objective is to find what are the max number of vehicles UPS could deploy for delivering packages to the chipotle locations in Dallas.

In [11]:
import pandas as pd
num_vehicles = [1,2,3,4,5]

results_df = pd.DataFrame(columns=['numberofvehicles','maxdistance'])

for veh_count in num_vehicles:
    result  = opt_vrp(veh_count)
    results_df = results_df.append({'numberofvehicles':veh_count,
                                     'maxdistance':result},
                                      ignore_index = True)
q,r: (1, 6, 13)
<ortools.constraint_solver.pywrapcp.RoutingIndexManager; proxy of <Swig Object of type 'operations_research::RoutingIndexManager *' at 0x7f1696730390> >
Objective: 7813360
Route for vehicle 0:
 0 ->  3 ->  2 ->  5 ->  4 ->  1 ->  6 ->  9 ->  7 ->  11 ->  12 ->  8 ->  10 -> 0
Distance of the route: 77360m

Maximum of the route distances: 77360m
q,r: (1, 6, 13)
<ortools.constraint_solver.pywrapcp.RoutingIndexManager; proxy of <Swig Object of type 'operations_research::RoutingIndexManager *' at 0x7f169673ab40> >
Objective: 4936795
Route for vehicle 0:
 0 ->  8 ->  7 ->  9 ->  6 ->  11 ->  12 ->  10 -> 0
Distance of the route: 42234m

Route for vehicle 1:
 0 ->  3 ->  5 ->  4 ->  1 ->  2 -> 0
Distance of the route: 48461m

Maximum of the route distances: 48461m
q,r: (1, 6, 13)
<ortools.constraint_solver.pywrapcp.RoutingIndexManager; proxy of <Swig Object of type 'operations_research::RoutingIndexManager *' at 0x7f169647cf30> >
Objective: 4341684
Route for vehicle 0:
 0 ->  8 ->  7 ->  9 ->  6 ->  11 ->  12 ->  10 -> 0
Distance of the route: 42234m

Route for vehicle 1:
 0 ->  2 ->  5 -> 0
Distance of the route: 40940m

Route for vehicle 2:
 0 ->  3 ->  4 ->  1 -> 0
Distance of the route: 35110m

Maximum of the route distances: 42234m
q,r: (1, 6, 13)
<ortools.constraint_solver.pywrapcp.RoutingIndexManager; proxy of <Swig Object of type 'operations_research::RoutingIndexManager *' at 0x7f169638c960> >
Objective: 4240090
Route for vehicle 0:
 0 ->  3 ->  4 ->  1 -> 0
Distance of the route: 35110m

Route for vehicle 1:
 0 ->  8 ->  7 ->  9 ->  12 -> 0
Distance of the route: 32494m

Route for vehicle 2:
 0 ->  2 ->  5 -> 0
Distance of the route: 40940m

Route for vehicle 3:
 0 ->  6 ->  11 ->  10 -> 0
Distance of the route: 37546m

Maximum of the route distances: 40940m
q,r: (1, 6, 13)
<ortools.constraint_solver.pywrapcp.RoutingIndexManager; proxy of <Swig Object of type 'operations_research::RoutingIndexManager *' at 0x7f1696285120> >
Objective: 4240090
Route for vehicle 0:
 0 ->  3 ->  4 ->  1 -> 0
Distance of the route: 35110m

Route for vehicle 1:
 0 ->  8 ->  7 ->  9 ->  12 -> 0
Distance of the route: 32494m

Route for vehicle 2:
 0 ->  2 ->  5 -> 0
Distance of the route: 40940m

Route for vehicle 3:
 0 -> 0
Distance of the route: 0m

Route for vehicle 4:
 0 ->  6 ->  11 ->  10 -> 0
Distance of the route: 37546m

Maximum of the route distances: 40940m

Plot the sensitivity using plotly :   Now plot a graph using plotly with x =number of vehicles y = max distance

In [13]:
import plotly.express as px
fig =px.scatter(results_df, x= "numberofvehicles", y="maxdistance",trendline ="lowess")
fig.show()

Conclusion

Though this project we have the following recommendations for UPS to deliver packages to the chipotle locations in Dallas, TX.
a. The max distance any vehicle should travel is 40940m
b. The max number of vehicles UPS should deploy is 4.